About one hundred years ago, it became possible to predict the thermal dynamical properties and phase transform conditions of realistic systems in gaseous, liquid, or solid states due to the birth of statistical physics. Up to the present, however, the predication is very limited because it is difficult to accurately compute the partition functions (PFs) of systems composed of more than several hundreds atoms. Since the PF can directly produce the state equation, the internal energy, the entropy, and the free energy, which represents the driving force of a system in nonequilibrium evolutions, precise calculations of the PF are highly desired today in nanoclusters synthesis,^[1,2] predicting protein folding,^[3–5] materials’ phase transformation,^[6] and designing new functional materials as crystal^[7,8] or bulk metal glass.^[9]

The attempts to exactly compute PFs and free energy have been continued in the past half a century. In 1950s, Zwanzig et al. put forward the earliest approach to calculate the free energy,^[10,11] which is able to obtain the free energy difference between two states. Nevertheless, the method is largely limited by the low accuracy and can only be cautiously applied in cases where disparity between the two states is small.^[12] Several methods were later developed by calculating the density of states, eg., histogram reweighting,^[13] and transition matrix.^[14] Among them, Wang–Landau (W–L) sampling formalism^[15] proved to be the most effective for its successful applications to small systems such as the Ising model^[15] and lattice protein model.^[16] However, the computational efficiency and accuracy deeply rely on the convergence criterion and the modification factor.^[17–20] In 2011, Hainam Do et al.^[21] greatly improved the W–L method and proposed a more efficient method, which on the same accuracy level, takes into account the amount of configurations about 2–3 orders less than the original W–L method.

In spite of the fact that considerable progress has been made in PF computations, it goes far beyond the current computer capacity to deal with macroscopic systems.^[20,22–24] The W–L method, for example, has to consider configurations even only concerning 64 water molecules.^[21] In fact, the upper limit treated by the W–L method has always been less than 300 particles and it could be hardly possible to extend the previous methods to larger systems. Though Hainam Do et al. improved the efficiency, there is still a huge amount of calculation when the system becomes larger. Furthermore, the accuracy of the previous methods is also a problem. Usually, it was tested by comparisons with experimental measurements of magnetization,^[15] specific heat capacity,^[17] or pressure.^[6] However, the theoretical results depend greatly on some artificial parameters in the empirical potential (EP), so even if the theoretical and experimental results are in good agreement, it cannot yet be affirmed that the method for calculating the PF is accurate because the EP may not be correct to describe the realistic systems. For example, two EPs were employed to calculate the pressure of a carbon dioxide system in a two-phase equilibrium state by the same method^[6] and produced significantly different results, one of which deviates obviously from the experimental measurements. Clearly, the deviation cannot deny the theoretical method for calculating the pressure. Relatively, molecular dynamics (MD) simulation^[25–29] should be a more reasonable way to examine the methods because the MD and the PF computation share the same EP. Actually, we have performed MD simulations of gaseous Lennard–Jones (L–J) atoms to obtain the internal energy and the isothermal expansion work, which were compared with the ones obtained from the PF using the same L–J potential.^[29] When the parameters of the L–J potential were adjusted to be so small that the atoms are nearly ideal gaseous atoms, the MD results are in very good agreement with the ones resulted from the PF of ideal gas obtained analytically, showing that the MD simulations are an excellent way for testing the PF. When the same MD procedure was performed on Ar atoms with L–J interaction potential, the internal energy and the expansion work at 200–300 K differ by about 20%–30% from the ones resulted from the PF obtained by the method of Hainam Do et al.^[21,30]

In this paper, we proposed a more rapid and accurate method to calculate the PF of canonical ensembles and applied it to systems containing up to 1000 Ar atoms in either liquid or dense gaseous states. The obtained internal energies and isothermal expansion works (or free energy) are in a good agreement with those from MD simulations. Lastly, we illustrated that the method works easily for macroscopic or even infinite systems.

For a canonical ensemble composed of systems at temperature T of volume V containing N indistinguishable particles of mass m with their coordinates and momenta denoted by and , respectively, the PF reads

(1)

(2)

(3)

Although departure of the configuration integral decreases the number of the integral variables in Eq. (1) by a half, we consider that it is this departure that leads to the difficulties existing in calculation of the PFs, because, for general potential function , it is too difficult to solve the 3N-fold integral neither by the analytical method nor via numerical calculation if the number of particles goes up to hundreds.

Different from the previous methods, we do not make the departure, but instead, write the integral (1) as

(4)

(5)

(6)

All the quantities in Eq. (6) can be easily obtained by MD and MC simulations. Specifically, every particle of the system is assigned a certain velocity initially and then is allowed to relax freely for hundreds of picoseconds to produce , , , and the potential range as functions of . For systems composed of up to 1000 Ar atoms confined in a cubic box, described by L–J interaction potential^[32] with and ε = 119.8 K, the MD simulations with an integral step of 1 fs need only to last about 400 ps to produce the quantities, and fluctuations of the total momentum, the total energy, and the total potential are indeed small (as shown in Fig. 1) for 500 and 864 Ar atoms in a density of or , which are in agreement with our above analysis, although in principle the kinetic energy can completely convert into the potential energy. In consideration of the fluctuation, all the sampled data such as shown in Fig. 1 were fitted by smooth functions of to produce , , , and , which were used in Eq. (6) to determine . For obtaining Q, one thousand configurations have been sampled and stored during the above MD simulations with given , and for each one of the configurations, an atom selected arbitrarily was randomly positioned (shooting) within the cubic box for ten thousands times ( ) while the coordinates of other atoms remained unchanged; and for each random shooting, if was valid, then the allowed number K increased by 1, and finally we obtained the possible space volume for the atom , and . In consideration of the fluctuation due to different configurations, we took the averaged value in Eq. (6) and obtained the state density, as shown in Fig. 2. Finally, the data were fitted to functions of (solid curves in Fig. 2) to be used to calculate the PFs, which produce the internal energy, the entropy, the free energy, and so on.

Fig. 1.

	Figure Option ViewDownloadNew Window
	Fig. 1. (color online) The variance range of the total momentum ( , ), the potential ( , ), , in arb. units as functions of for 500 Ar atoms in the density of (a) or (b) , and for 864 Ar atoms in the density of (c) or (d) . Here the total potential is shifted and is amplified for clear display.

Fig. 2.

	Figure Option ViewDownloadNew Window
	Fig. 2. (color online) The calculated state density (dots) of 500, 864, and 1000 Ar atoms in a density of (gaseous state, denoted by G) or (liquid state, denoted by L) and the fitted functions (solid curves). Note that the independent variable for G1000 and L864 is instead of .

In order to test the accuracy of the PFs, the same L–J potential was applied in MD simulations of the system contacted with a heat reservoir at temperature T to produce the internal energy and isothermal expansion work at temperature T to compare with the results of the PFs. To obtain the internal energy, all atoms of the system were assigned to velocities of Maxwell’s distribution at the temperature T every 100 fs in the early period of 1 ps of the MD simulation with time step of 1 fs, and then the system was allowed to evolve for 400 ps to produce the mean value of the internal energy, during which ten atoms randomly selected every 0.1 ps were assigned to velocities of the Maxwell’s distribution to simulate effects of the heat reservoir.^[25–29] For obtaining the expansion work from the MD simulations, similar MD simulations were performed for different volumes of the system. Specifically, one side length of the cubic box was enlarged by a step of 0.1 Å, until the length increases by 2 Å, and for each step, the MD simulation lasted 400 ps to produce the pressure F calculated by^[32]

(7)

The internal energy and the expansion work at different temperature T from MD simulations for the systems composed of 500, 864, 1000 Ar atoms in a density of (gaseous state) or (liquid state) are shown in Fig. 3, where the results from the PFs of our method are also presented. The comparisons show that the differences between the MD and PFs in both the internal energy and the expansion work are no more than 3% for lower temperature, and smaller than 1% for high temperature. It is notable that the gas density of employed in our calculation is about three hundred times of the atmospheric density ( ) under standard conditions and the concerned pressure reaches up to nine hundreds times of one atmospheric pressure, while the density of corresponds to a liquid state and the pressure concerned in our calculation is as high as Pa, which is nearly ten thousands times of one atmospheric pressure.

Fig. 3.

	Figure Option ViewDownloadNew Window
	Fig. 3. (color online) (a) The internal energy and (b) the expansion work at different temperature T from MD simulations (dot) and our method (solid line) for the systems composed of 500, 864, 1000 Ar atoms in a density of (gaseous state, denoted by G) or (liquid state, denoted by L). The internal energy for the liquid state was shifted by 60 eV for clear displaying.

For comparison, we also calculated the PFs of a system composed of 500 Ar atoms in a density of by MHD to produce the internal energy and the expansion work at different temperatures, which are presented in Table 1, Table 2, and Fig. 4 to be compared with the results of the MD. The relative deviation in both the internal energy and the expansion work is obviously larger than the one from our method, especially at lower temperatures (200–300 K). As an example, the internal energy and the expansion work at 200 K from the MHD are 2.6 eV and , respectively, while the results of the MD are 3.7 eV and , resulting in relative deviations as large as 30% and 20%, respectively, which are about ten times larger than the ones (∼3%) of our method. It should be stated that the accuracy of the MHD depends on the amount of the atomic configurations sampled in the calculations using the MHD, and 30000 configurations, which was suggested in Ref. [21], were sampled in each step to ensure that further increasing the number of configurations does not change the results significantly.^[30] As for the computing speed, the MHD is about 10 times slower than our method, that is why we did not apply the MHD to calculate the PFs of larger systems.

Table 1.

Table 1.

Table 1. The internal energy (in eV) at different temperature T (in K) from MD simulations, our method, and MHD for the system composed of 500 Ar atoms in a density of . . Temperature MD Ours MHD 200 3.7 3.8 2.6 300 11.7 11.2 10.0 400 18.5 18.4 17.2 500 25.5 25.7 24.2 600 32.1 32.2 31.2 700 39.8 40.2 38.2 800 46.0 46.2 45.2

Table 1. The internal energy (in eV) at different temperature T (in K) from MD simulations, our method, and MHD for the system composed of 500 Ar atoms in a density of . .

Table 2.

Table 2.

Table 2. The expansion work (in ) at different temperature T (in K) from MD simulations, our method, and MHD for the system composed of 500 Ar atoms in a density of . . Temperature MD Ours MHD 200 1.77 1.70 1.38 300 3.26 3.20 3.07 400 4.80 4.73 4.72 500 6.20 6.31 6.32 600 7.82 7.85 7.92 700 9.07 9.27 9.46 800 10.8 10.6 11.0

Table 2. The expansion work (in ) at different temperature T (in K) from MD simulations, our method, and MHD for the system composed of 500 Ar atoms in a density of . .

Fig. 4.

	Figure Option ViewDownloadNew Window
	Fig. 4. (color online) (a) The internal energy and (b) the expansion work at different temperature T from MD simulations (dot), our method (solid line), and MHD (dash line) for the systems composed of 500 Ar atoms in a density of . The relative deviations of MHD (red dot) and our method (black dot) are shown in the inset.

It should be noted that the state density increases with the total energy E by , while the Boltzmann factor decreases exponentially with E for given T, resulting in a sharp distribution of the total energy ( ) around the most probable energy , corresponding to the temperature T. That is, the state density at E distant from by has little contribution to the PF for the given temperature T. So, for studying a system at given temperature via calculating the PF, we only need to calculate the state density around , which can be easily implemented by our method to further save the computation task. This makes our method much faster than MHD, which has to calculate the state density of the potential energy ranged from the highest to the lowest. Thus, our method can easily produce the PFs of much larger systems composed of millions of atoms that are large enough to ensure that the internal energy and the free energy increase linearly with the number of the atoms in the systems. As shown in Fig. 5, for the Ar systems, the energy ε and the free energy f of a single atom are indeed constants if the systems contain more than 500 atoms. That is, our method is applicable to macroscopic systems to produce the internal energy and the free energy by and .

Fig. 5.

	Figure Option ViewDownloadNew Window
	Fig. 5. (color online) (a) The energy ε and (b) the free energy f per atom for systems composed of 500, 680, 864, and 1000 Ar atoms in density of for different temperatures.

Recognizing the flaw of departing the configurational integral from the PFs integral, we start from calculation of the state density of the total energy rather than the total potential energy to calculate the PFs, showing that it is much more accurate and efficient than the previous methods and may be applicable to macroscopic systems. Obviously, our method would find its vast application areas immediately.